
EE.value <- function(data, cnt, alpha.start,beta, gamma,
                                Hessian=FALSE,y.hat){
  p.alpha   = length(alpha.true);   p.beta = length(beta.true)
  startpars <- alpha.hat <- c(alpha.start)
  # pars <- startpars <-alpha.hat <- alpha.true # DEBUG
  # beta <- beta.true; gamma <- gamma.true #DEBUG
  l0 = data[,c("l0.1","l0.2")]
  a0 = data[,"a0"];  l1 = data[,"l1"];  a1 = data[,"a1"];   y = data[,"y"]
  nb = nrow(data)  
  
  alpha.hat <- matrix(alpha.hat[1:p.alpha],ncol=2);
  theta.hat = exp(l0 %*% t(alpha.hat))
  phi123 = exp(l0 %*% t(beta[1:3,]))
  phi45  = expit(l0 %*% t(beta[4:5,]))
  p.l1 = sapply(1:nb, function(i) phi45[i,2 - a0[i]])
  # Nuisance functions
  p.a1 = c(expit(cbind(l1,a0,l0) %*% gamma["a1",]))
  p.a0 = c(expit(l0 %*% gamma["a0",1:2]))
  k = 2 + 2*(1-a0) + (1-l1)
  
  theta.k.hat <- sapply(1:nb, function(i) theta.hat[i,k[i]])
  theta.k.inv.hat <- 1/theta.k.hat
  theta.k.inv.a1.hat <- (theta.k.inv.hat)^a1
  
  theta..hat <- sapply(1:nb, function(i) theta.hat[i,1])
  theta..inv.hat <- 1/theta..hat
  theta..inv.a0.hat <- (theta..inv.hat)^a0
  
  dr.objective <- function(pars){
    # pars represents alpha here. beta is beta, gamma is gamma.
    # base.covs <- func.DataGen(c(pars, c(beta), c(gamma)), data,cnt)
    # MyAttach(base.covs)
    pars <- matrix(pars[1:p.alpha],ncol=2);
    theta  = exp(l0 %*% t(pars)); 
    
    theta.k <- sapply(1:nb, function(i) theta[i,k[i]])
    theta.k.inv <- 1/theta.k
    theta.k.inv.a1 <- (theta.k.inv)^a1
    
    # for (i in 1:nb){theta.k[i] <- theta[i,][k[i]]}
    
    ### m = 1
    u1 <- y * theta.k.inv.a1 
    if (ee.assume.true.alpha){
      E.u1 <- g(l1, a0, l0, y.hat)
    }else{
      E.u1 <- g(l1, a0, l0, p.y.cond, p.a1, 
                theta.k.inv=theta.k.inv.hat, gamma, nb,y.hat)
      # print(1) #DEBUG
    }
    
    # print(u1 - E.u1)
    
    # d1.row.k <- -l0 * a1 * p.y * theta.k.inv.a1
    if (ee.assume.true.alpha){
      d1.row.k <- -l0 * a1 * y.hat.func(l0,a0,l1,a1=vec.0,y.hat)
      E.d1.row.k <- -l0 * g1(l1, a0, l0, y.hat,gamma)
    }else{
      d1.row.k <- -l0 * a1 * y.hat.func(l0,a0,l1,a1,y.hat) * theta.k.inv.a1.hat
      # E.d1.row.k <- -l0 * p.a1 * E.u1
      # E.d1.row.k <- -l0 * g1(l1, a0, l0, p.y.cond, p.a1, theta.k.inv, gamma, nb)
      E.d1.row.k <- -l0 * g1(l1, a0, l0, theta.k.inv=theta.k.inv.hat, gamma, nb,y.hat)
      # print(2) # DEBUG
    }
    d1 <- array(rep(0,5*2*nb),dim=c(nb,5, 2))
    for (i in 1:nb) {
      d1[i,k[i],] <- d1.row.k[i,]
    }
    E.d1 <- array(rep(0,5*2*nb),dim=c(nb,5, 2))
    for (i in 1:nb) {
      E.d1[i,k[i],] <- E.d1.row.k[i,]
    }
    
    # test.d1 <- 0 #DEBUG
    # for (i in 1:nb){test.d1 <- test.d1+((d1 - E.d1))[i,,] * cnt[i]}
    # 
    # test.u1 <- 0 #DEBUG
    # for (i in 1:nb){test.u1 <- test.u1+((d1 - E.d1) * (u1 - E.u1))[i,,] * cnt[i]}
    # print(test.u1)
    # j <<- j+1
    # vals[[j]] <<- test.u1
    
    ### m = 0
    theta. <- theta[,1] #OLD: sapply(1:nb, function(i) theta[i,1])
    theta..inv <- 1/theta.
    theta..inv.a0 <- (theta..inv)^a0
    
    u0 <- y * theta.k.inv.a1 * theta..inv.a0
    p.a0 <- p.a0.func(l0,gamma)
    E.u0 <- p.a0 * h(a0=vec.1, l0, theta..inv.a0= theta..inv.a0.hat,
                     p.y.cond, p.a1, theta.k.inv=theta.k.inv.hat, 
                     gamma, nb,
                     beta, y.hat,l1.hat) +
      (1-p.a0) * h(a0=vec.0, l0, theta..inv.a0= theta..inv.a0.hat,
                   p.y.cond, p.a1, theta.k.inv=theta.k.inv.hat, gamma, nb,
                   beta, y.hat,l1.hat)
    # test.u0 <- 0 #DEBUG
    # for (i in 1:nb){test.u0<- test.u0+((u0 - E.u0))[i] * cnt[i] }
    # print(test.u0) #DEBUG
    # j <<- j+1
    # vals[[j]] <<- test.u1
    
    
    d0.row.1 <- f.(a0,l0, theta..inv.a0= theta..inv.a0.hat, 
                   p.y.cond, p.a1, theta.k.inv=theta.k.inv.hat
                   , gamma, nb,
                   beta,y.hat,l1.hat)
    # theta..inv.a0.hat <- exp(exp(l0 %*% t(alpha.hat)))-# DEBUG
    d0.row.k <- fk(a0,l0, theta..inv= theta..inv.hat,
                   p.y.cond, p.a1, theta.k.inv=theta.k.inv.hat, gamma, nb,
                   beta,y.hat,l1.hat) # This is a matrix of dim = (4, nb)
    
    # p.a0 / theta..inv.hat + (1-p.a0) #DEBUG
    # E.d0.row.1: Two should be equivalent # Can be sanity check
    E.d0.row.1.test <- p.a0 * f.(a0=vec.1,l0, theta..inv.a0,
                                 p.y.cond, p.a1, theta.k.inv, gamma, nb,
                                 beta, y.hat,l1.hat) +
      (1- p.a0) * f.(a0=vec.0,l0, theta..inv.a0,
                     p.y.cond, p.a1, theta.k.inv, gamma, nb,
                     beta, y.hat,l1.hat)
    E.d0.row.1 <- -l0 * p.a0 * h(a0=a0, l0, 
                                 theta..inv.a0= theta..inv.a0.hat,
                                 p.y.cond, p.a1, theta.k.inv=theta.k.inv.hat, 
                                 gamma, nb,
                                 beta, y.hat,l1.hat)
    
    # E.d0.row.k <- p.a0 * fk(a0=vec.1,l0, theta..inv= theta..inv.hat,
    #                         p.y.cond, p.a1, theta.k.inv=theta.k.inv.hat, gamma, nb,
    #                         beta, y.hat,l1.hat) +
    #   (1- p.a0) * fk(a0=vec.0,l0, theta..inv= theta..inv.hat,
    #                  p.y.cond, p.a1, theta.k.inv=theta.k.inv.hat, gamma, nb,
    #                  beta, y.hat,l1.hat)
    E.d0.row.k <- E.fk(p.a0=p.a0,l0, theta..inv= theta..inv.hat,
                       p.y.cond, p.a1, theta.k.inv=theta.k.inv.hat, gamma, nb,
                       beta,y.hat,l1.hat) # This is a matrix of dim = (4, nb)
    
    
    d0 <- array(rep(0,5*2*nb),dim=c(nb,5, 2))
    for (i in 1:nb) {
      d0[i,1,] <- d0.row.1[i,]
      # d0[i,k[i],] <- d0.row.k[i,]
      # d0[i,2:5,] <- d0.row.k[i,]
    }
    # change 27Jul2021
    d0[,2:5,] <- d0.row.k
    
    E.d0 <- array(rep(0,5*2*nb),dim=c(nb,5, 2))
    for (i in 1:nb) {
      E.d0[i,1,] <- E.d0.row.1[i,]
      # E.d0[i,k[i],] <- E.d0.row.k[i,]
      # E.d0[i,2:5,] <- E.d0.row.k[i,]
    }
    # change 27Jul2021
    E.d0[,2:5,] <- E.d0.row.k
    
    tmp <-  array(rep(0,5*2),dim=c(5, 2))
    tmp.stage0 <- tmp.stage1 <-  array(rep(0,5*2),dim=c(5, 2)) # DEBUG
    tmp.d1 <- tmp.d0 <-  array(rep(0,5*2),dim=c(5, 2)) # DEBUG
    tmp.u1 <- tmp.u0 <- 0
    # tmp2 <- matrix(rep(0,32*2),ncol=2) #  DEBUS
    for (i in 1:nb){
      total <- (d1 - E.d1) * (u1 - E.u1) + (d0 - E.d0) * (u0 - E.u0)
      # total <- (d1 - E.d1) * (u1 - E.u1)
      # tmp <- tmp + (total[i,,])^2 * cnt[i]
      tmp <- tmp + total[i,,] * cnt[i]
      tmp.stage1 <- tmp.stage1 + ((d1 - E.d1) * (u1 - E.u1))[i,,] * cnt[i]
      tmp.d1 <- tmp.d1 + ((d1 - E.d1))[i,,] * cnt[i]
      tmp.u1 <- tmp.u1 + ((u1 - E.u1))[i] * cnt[i] 
      tmp.stage0 <- tmp.stage0 + ((d0 - E.d0) * (u0 - E.u0))[i,,] * cnt[i]
      tmp.d0 <- tmp.d0 + ((d0 - E.d0))[i,,] * cnt[i]
      tmp.u0 <- tmp.u0 + ((u0 - E.u0))[i] * cnt[i] 
      # tmp2 <- ((d1 - E.d1) * (u1 - E.u1))[i,,] * cnt[i]
      # tmp3 <- tmp3 + (u0 - E.u0)[i]
      # tmp2 <- tmp2 + (d0.row.1  - E.d0.row.1 )*cnt[i]
    }
    
    # return(sum(tmp1))
    # return(sum(tmp))
    return(list(tmp=tmp, tmp.stage1=tmp.stage1,
                tmp.d1=tmp.d1, tmp.u1=tmp.u1,
                tmp.d0=tmp.d0, tmp.u0=tmp.u0))
  }
  
  return(dr.objective(c(startpars)))
}